function welfare = Ramsey_plus_affine_without_epsilon_grid(pr,w,pai,pw,G, T_grid, tau_grid)

% Parameters and Grids 

% Parameters
% pr.sig = 1/0.5;              % inverse Frisch
% pr.tau = 0.181 ;               % benchmark progressivity
% h_0_target = 0.3;           % target for labor supply without transfers  % do not need now
% phi = 1;                    %  from HT clibration 
% G =   0.181248550503625;    %  from HT clibration 
% rhoa =1;                    % persistence idiosyncratic shock

% % Targets
% v_log_con = 0.332*(1-0.296);                            % variance log consumption data
% v_log_inc = 0.4117;                                     % variance log income data
% vw   = (1-rhoa^2)*v_log_con*(1/(1-pr.tau)^2);              % variance uninsurable shock
%    

T_n=length(T_grid);
tau_n=length(tau_grid);

% Stack all combinations of taxes and transfers
fiscal = gridmake(tau_grid,T_grid);




%%%%%%%%%% From HT %%%%%%%%% 

% % Uninsurable shocks
% Discertizing as normal
% nodes_alpha_n = 100;
% alpha_mean = -vw/(2*(1-rhoa^2));
% alpha_var = vw/(1-rhoa^2);
% [nodes_alpha,weights_alpha] = qnwnorm(nodes_alpha_n,alpha_mean,alpha_var);

% %% Preference
% pr.gam = 1;                         % risk aversion
% pr.sig = pr.sig;                         % Frisch elasticity

% %% Government 
% pr.g     = 0.188;                   % fraction of Government Expenditure
% pr.theta = 0;                       % taste for redistribution

% %% Wage Distribution: Setup
% pr.na     = 500;                  % number of grid points for alpha shocks
% nodes_alpha_n=pr.na 
% pr.wmin   = 0.05;                   % minimum wage  
% pr.wmax   = 74;                     % maximum wage
% pr.pareto = 2.2;                    % Pareto tail parameter 
% pr.vy     = 0.4117;                 % normal variance of log earnings
% pr.vc     = 0.332*(1-0.296);        % variance of log consumption
% 
% pr.normal_alpha= pr.vc   / (1-pr.tau)^2 - 1/ pr.pareto^2;
% veps=  (pr.vy- pr.normal_alpha)*(pr.sig/(pr.sig+1))^2;
%  
% % Discretization alpha
% pr.va     = pr.vc/(1-pr.tau)^2 ;    %  -1/pr.pareto^2                               % variance of uninsurable shock
% pr.ve     = (pr.vy-(pr.va-1/pr.pareto^2))*((1+pr.sig)/pr.sig)^-2; % variance of insurable shocks
% temp.ans  = fsolve(@(x) alphadist(x,pr),[log((pr.pareto-1)/pr.pareto)-(pr.va-1/pr.pareto^2)/2 (pr.va-1/pr.pareto^2)^(1/2)],optimset('TolX',1e-12,'Display','off'));  % alpha distribution
% alphadist(temp.ans,pr); clc

%load w                            % exp(alpha)
nodes_alpha=log(w);
weights_alpha=pai;% probablity measure

%load pai
%load pdf      

% Insurable shocks  
%veps=(pr.vy-(pr.va-1/pr.pareto^2))*((1+pr.sig)/pr.sig)^-2;
% We do not need, note that I calculated labor
% disutility without epsilon
nodes_eps_n = 1;
% eps_mean = -veps/2;
% eps_var = veps;

%[nodes_eps,weights_eps] = qnwnorm(nodes_eps_n,eps_mean,eps_var);
nodes_alpha_n=length(w);

% Stack all combinations of insurable and uninsurable shocks

%% Labor Supply, Consumption, Government Budget

tol = 1e-8;
max_iter = 5000;
lambda_grid = zeros(T_n*tau_n,1);
ymean = zeros(T_n*tau_n,1);
c_mat = zeros(nodes_alpha_n,T_n*tau_n);
H_mat = zeros(nodes_alpha_n,T_n*tau_n);
labor_disutility = zeros(nodes_alpha_n,T_n*tau_n);
util_mat=zeros(nodes_alpha_n,T_n*tau_n);
% n_mat = zeros(nodes_alpha_n,nodes_eps_n,T_n*tau_n);
% y_mat = zeros(nodes_alpha_n,nodes_eps_n,T_n*tau_n);
% yhat_mat = zeros(nodes_alpha_n,nodes_eps_n,T_n*tau_n);
% 
% y_mat_family = zeros(nodes_alpha_n,T_n*tau_n);
% c_mat_family = zeros(nodes_alpha_n,T_n*tau_n);


omega=exp((1+pr.sig)/pr.sig^2 * pr.ve/2); % note that this is an expectation of E exp((1+pr.sig)/pr.sig epsilon)

for i = 1:T_n*tau_n
    lambda_min = 0.2;
    lambda_max = 1.2;

    disp(['Fiscal system indexed ', num2str(i)])

    
    
    for iter_fiscal = 1:max_iter
        
        
       if iter_fiscal == max_iter
            disp('Maximum number of iterations reached for lambda')
            disp(['Fiscal system indexed ', num2str(i)])
        end 
        
        lambda = 0.5*(lambda_min+lambda_max);
        
            
        H_min = 0 + 0*nodes_alpha;
        H_max = 100 + 0*nodes_alpha;

        for iter_alpha = 1:max_iter
                        
            H = 0.5*(H_min+H_max);
            nominator = lambda^(1/pr.sig)*(1-fiscal(i,1))^(1/pr.sig)*exp(nodes_alpha*(1-fiscal(i,1))/pr.sig) .*H.^(-fiscal(i,1)/pr.sig);
            if pr.gam ==1
                denominator=(lambda * exp(nodes_alpha*(1-fiscal(i,1))) .*H.^(1-fiscal(i,1)) + fiscal(i,2)).^(1/pr.sig);
            else
                denominator=(lambda * exp(nodes_alpha*(1-fiscal(i,1))) .*H.^(1-fiscal(i,1)) + fiscal(i,2)).^(pr.gam/sigma);
            end
            err_H = H - nominator .* omega./denominator;
%    
            if max(abs(err_H))< tol
                H_mat(:,i) = H;
                c_mat(:,i) = lambda*exp(nodes_alpha*(1-fiscal(i,1))).*H.^(1-fiscal(i,1))+fiscal(i,2)  ;
                labor_disutility(:,i) =(nominator./denominator).^(1+pr.sig)*omega/(1+pr.sig);  
                y_mat(:,i) = exp(nodes_alpha).*H;
                yhat_mat(:,i) = lambda*y_mat(:,i).^(1-fiscal(i,1));
%                     y_mat_family(j,i) = exp(nodes_alpha(j))*H;
%                     c_mat_family(j,i) = lambda*exp(nodes_alpha(j)*(1-fiscal(i,1)))*H^(1-fiscal(i,1))+fiscal(i,2);
                
                break
            else
                
                H_max =H.*(err_H>0) + H_max.*(err_H<=0);
                H_min =H_min.*(err_H>0) + H.*(err_H<=0);
            end
                
    
        end         
%             
%             for k = 1:nodes_eps_n
%               n_mat(j,k,i) =(nominator/denominator)*exp(nodes_eps(k)*1/pr.sig) ;
%               y_mat_ind(j,k,i) =(nominator/denominator)*exp(nodes_eps(k)*1/pr.sig)*exp(nodes_alpha(j)+nodes_eps(k))  ;       
%             end
        
           
    Y = sum(pai.*y_mat(:,i));
    Yhat = sum(pai.*yhat_mat(:,i));
    ymean(i) = Y;
         
        err_fiscal = G  +fiscal(i,2)- (Y-Yhat) ;% + fiscal(i,2)
   
        if abs(err_fiscal) < tol
            lambda_grid(i) = lambda;
            break
        elseif err_fiscal < 0
            lambda_min = lambda;
        else
            lambda_max = lambda;
        end
     end
   
end

% analytically  (1+pr.sig) /pr.sig* veps+pr.va+1/pr.pareto^2
if min(min(min(c_mat))) < 0
    disp('Consumption negative')
end
% if min(min(min(n_mat))) < 0
%     disp('Labor negative')
% end

%% Welfare
if pr.gam==1
util_mat  =log(c_mat) -labor_disutility;
else
util_mat  =c_mat.^(1-pr.gam)/(1-pr.gam) -labor_disutility;
end

W = zeros(T_n*tau_n,1);
for i = 1:T_n*tau_n
    W(i) = sum(pai.*util_mat(:,i));
end
                     
[value,loc] = max(W);
welfare=value;

tau = fiscal(loc, 1);
Tr   = fiscal(loc, 2);
lambda=lambda_grid(loc);

c=c_mat(:,loc);
y=y_mat(:,loc);
save lambda lambda
save tau tau
save Tr Tr
save c c
save y y 

